#!/usr/bin/perl
# 2012/10/17 Calculate distance from nearest restriction sites

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [bowtie aligned file]\n";
}

my %opt;
getopts("i:", \%opt);
my $FILE_bowtie = $opt{i};

#---------------------------------------
# read restriction file
#---------------------------------------
my %num2pos;	# fragment numberから具体的な位置を示す
my %cate2num;	# categoryからそこに属するfragment numberの一群を示す

#number	chr	position	length_before	length_after
#100001	chr1	1004	1004	230
#100002	chr1	1234	230	1151
#100003	chr1	2385	1151	822
#100004	chr1	3207	822	30
my $FILE_enzyme = '/home/hideki/Project/2012-10-16_DamID_KD_swi6_cut14/out/DpnI_sites.txt';
my $fh_en = IO::File->new($FILE_enzyme) or die "cannot open $FILE_enzyme: $!";
$fh_en->getline();
while($_ = $fh_en->getline()){
	s/\r?\n//;
	my ($number, $chr, $position, $lenB, $lenA) = split /\t/;
	my $cate = int($position / 1000) * 1000;
	$num2pos{$number} = $position;
	push @{$cate2num{$cate}}, $number;
}
$fh_en->close();

#---------------------------------------
# read target file
#---------------------------------------
# 断片の端っこがlocationの位置
#ILLUMINA-CDA390:98:FC:3:1:17901:1020    +       chr3    2099840 1
#ILLUMINA-CDA390:98:FC:3:1:1504:1020     -       chr3    11248   1
#ILLUMINA-CDA390:98:FC:3:1:7858:1020     +       chr3    10890   1
#ILLUMINA-CDA390:98:FC:3:1:15874:1020    -       chr1    4055757 1
#ILLUMINA-CDA390:98:FC:3:1:16453:1021    +       chr3    7572    1
#ILLUMINA-CDA390:98:FC:3:1:6192:1021     -       chr3    2442366 1
#ILLUMINA-CDA390:98:FC:3:1:17211:1021    -       chr2    3012381 1
#ILLUMINA-CDA390:98:FC:3:1:18089:1021    -       chr3    8530    1
#ILLUMINA-CDA390:98:FC:3:1:13771:1021    -       chr3    17556   1


sub searchList{
	my ($cate, $ref_cate, $ref_total) = @_;
	if($$ref_total != 0){
		return;
	}
	if(exists $cate2num{$cate}){
		$ref_cate->{$cate} = 1;
		$$ref_total++;
	}else{
		my $cate_plus = $cate + 1000;
		&searchList($cate_plus, $ref_cate, $ref_total);
		my $cate_minus = $cate - 1000;
		&searchList($cate_minus, $ref_cate, $ref_total);
	}
}


my $fh_bowtie = IO::File->new($FILE_bowtie) or die "cannot open $FILE_bowtie: $!";
while($_ = $fh_bowtie->getline()){
	s/\r?\n//;
	my ($id, $direction, $chr, $position, $tmp) = split /\t/;
	my $cate = int($position / 1000) * 1000;

	my %list;
	my $total = 0;
	&searchList($cate, \%list, \$total);

	my @list;
	foreach my $c(keys %list){
		push @list, @{$cate2num{$c}};
	}

	my @sorted = sort {$a <=> $b } @list;
	my $MAX_plus1 = $sorted[$#sorted] + 1;
	my $min_minus1 = $sorted[0] - 1;

	my $min_distance = 9999999;
	foreach my $f(@sorted, $MAX_plus1, $min_minus1){
		if(exists $num2pos{$f}){
			my $restriction_site = $num2pos{$f};
			my $distance = 0;
			if($direction eq '+'){
				$distance = $position - $restriction_site;
			}elsif($direction eq '-'){
				$distance = $restriction_site - $position;
			}

			if($distance < 0){
				next;
			}

			if($distance < $min_distance){
				$min_distance = $distance;
			}
		}
		next;
	}
	print "$min_distance\n";
}
$fh_bowtie->close();

